#!/usr/bin/env python

"""yh-most-distant-sequence tries to find the 'most distant' sequence
according to a set of sequences, by alignment optimization.
"""

#This software is licensed under the MIT license:
#
# Copyright (c) 2013 ing. Youri Hoogstrate
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following
# conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
# WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.
#
#
# Documentation as defined and generated by:
# <http://epydoc.sourceforge.net/manual-fields.html#fields-synonyms>



__version_info__ = ["1","0","1","beta"]
__version__ = "-".join([".".join(__version_info__[:3])]+__version_info__[3:])
__author__ = "ing. Youri Hoogstrate"
__homepage__ = "https://code.google.com/p/yh-most-distant-sequence/"
__license__ = "MIT"
__organization__ = "Leiden University Medical Center: BioSemantics Group"



import random,math,datetime,time,sys



# Define constants:
INFINITY = 9999999999999



class SubstitutionMatrix:
	"""Controls substitution scores which can be used by the DistanceAlign class."""
	def __init__(self):
		"""Initiates the class by creating empty variables."""
		
		self.disabled = True
		self.scoring_matrix = {}
	
	def set_equal_scoring_probabilities(self,match = 1,mismatch_penalty = 0,gap_penalty = -0.5,alphabet = "ARNDCEQGHILKMFPSTWYV-"+"BZJXUO"):
		"""Sets the scoring of the substitution matrix to evolutionary neutral. The following scores are used by default:
			- Match penalty (any amino acid to only itself) = 1
			- Mismatch penalty (any amino acid to any other amino acid than itself) = 0
			- Gap penalty = -1.1 (This strictness is to avoid big gaps)
		
		---
		@param match: The score that should be obtained for a match
		@param mismatch_penalty: The mismatch penalty score; matching two unequal amino acids
		@param gap_penalty: The gap penalty for matching an amino acid with "-"
		@param alphabet: The alphabet of amino acids that the substitution matrix has.
		"""
		
		# Clean-up possible previous scoring matrix:
		del(self.scoring_matrix)
		self.scoring_matrix = {}
		
		# Fill the matrix:
		for char1 in alphabet:
			self.scoring_matrix[char1] = {}
			for char2 in alphabet:
				if(char1 == "-" or char2 == "-"):
					self.scoring_matrix[char1][char2] = gap_penalty
				elif(char1 == char2):
					self.scoring_matrix[char1][char2] = match
				else:
					self.scoring_matrix[char1][char2] = mismatch_penalty
		
		# Also remember the penalty values:
		self.mismatch_penalty = mismatch_penalty
		self.gap_penalty = gap_penalty
		
		self.disabled = False
	
	def set_blosum62(self,mismatch_penalty=0,gap_penalty=-0.5,alphabet = "ARNDCEQGHILKMFPSTWYV-"+"BZJXUO"):
		"""Sets the (normalized_ BLOSUM62 values as substitution scores. The original BLOSUM62 matrix is as follows:
		
		*   A   R   N   D   C   Q   E   G   H   I   L   K   M   F   P   S   T   W   Y   V \n
		A   4  -1  -2  -2   0  -1  -1   0  -2  -1  -1  -1  -1  -2  -1   1   0  -3  -2   0 \n
		R  -1   5   0  -2  -3   1   0  -2   0  -3  -2   2  -1  -3  -2  -1  -1  -3  -2  -3 \n
		N  -2   0   6   1  -3   0   0   0   1  -3  -3   0  -2  -3  -2   1   0  -4  -2  -3 \n
		D  -2  -2   1   6  -3   0   2  -1  -1  -3  -4  -1  -3  -3  -1   0  -1  -4  -3  -3 \n
		C   0  -3  -3  -3   9  -3  -4  -3  -3  -1  -1  -3  -1  -2  -3  -1  -1  -2  -2  -1 \n
		Q  -1   1   0   0  -3   5   2  -2   0  -3  -2   1   0  -3  -1   0  -1  -2  -1  -2 \n
		E  -1   0   0   2  -4   2   5  -2   0  -3  -3   1  -2  -3  -1   0  -1  -3  -2  -2 \n
		G   0  -2   0  -1  -3  -2  -2   6  -2  -4  -4  -2  -3  -3  -2   0  -2  -2  -3  -3 \n
		H  -2   0   1  -1  -3   0   0  -2   8  -3  -3  -1  -2  -1  -2  -1  -2  -2   2  -3 \n
		I  -1  -3  -3  -3  -1  -3  -3  -4  -3   4   2  -3   1   0  -3  -2  -1  -3  -1   3 \n
		L  -1  -2  -3  -4  -1  -2  -3  -4  -3   2   4  -2   2   0  -3  -2  -1  -2  -1   1 \n
		K  -1   2   0  -1  -3   1   1  -2  -1  -3  -2   5  -1  -3  -1   0  -1  -3  -2  -2 \n
		M  -1  -1  -2  -3  -1   0  -2  -3  -2   1   2  -1   5   0  -2  -1  -1  -1  -1   1 \n
		F  -2  -3  -3  -3  -2  -3  -3  -3  -1   0   0  -3   0   6  -4  -2  -2   1   3  -1 \n
		P  -1  -2  -2  -1  -3  -1  -1  -2  -2  -3  -3  -1  -2  -4   7  -1  -1  -4  -3  -2 \n
		S   1  -1   1   0  -1   0   0   0  -1  -2  -2   0  -1  -2  -1   4   1  -3  -2  -2 \n
		T   0  -1   0  -1  -1  -1  -1  -2  -2  -1  -1  -1  -1  -2  -1   1   5  -2  -2   0 \n
		W  -3  -3  -4  -4  -2  -2  -3  -2  -2  -3  -2  -3  -1   1  -4  -3  -2  11   2  -3 \n
		Y  -2  -2  -2  -3  -2  -1  -2  -3   2  -1  -1  -2  -1   3  -3  -2  -2   2   7  -1 \n
		V   0  -3  -3  -3  -1  -2  -2  -3  -3   3   1  -2   1  -1  -2  -2   0  -3  -1   4 \n
		
		The normalization finds the lowest value of the matrix (-4) and sets this is the value of 0 (no similarity at all).
		Next, the center the originial matrix, 0, is scaled to 1. The top, 11, is scaled to 3.75: new score = (score/-lowest score)+1
		
		---
		@param mismatch_penalty: The score for an amino acid if it isn't found in the BLOSUM62 matrix. Should be avoided because the substitution matrix is incomplete.
		@param gap_penalty: The score for extending or opening a gap.
		@param alphabet: The alplhabet of amino acids that are found in the biological sequences library.
		"""
		
		# Use lower triangle for filling.
		data = {'A':{'A':4,'C':0,'E':-1,'D':-2,'G':0,'F':-2,'I':-1,'H':-2,'K':-1,'M':-1,'L':-1,'N':-2,'Q':-1,'P':-1,'S':1,'R':-1,'T':0,'W':-3,'V':0,'Y':-2},'C':{'C':9,'E':-4,'D':-3,'G':-3,'F':-2,'I':-1,'H':-3,'K':-3,'M':-1,'L':-1,'N':-3,'Q':-3,'P':-3,'S':-1,'R':-3,'T':-1,'W':-2,'V':-1,'Y':-2},'E':{'E':5,'G':-2,'F':-3,'I':-3,'H':0,'K':1,'M':-2,'L':-3,'N':0,'Q':2,'P':-1,'S':0,'R':0,'T':-1,'W':-3,'V':-2,'Y':-2},'D':{'E':2,'D':6,'G':-1,'F':-3,'I':-3,'H':-1,'K':-1,'M':-3,'L':-4,'N':1,'Q':0,'P':-1,'S':0,'R':-2,'T':-1,'W':-4,'V':-3,'Y':-3},'G':{'G':6,'I':-4,'H':-2,'K':-2,'M':-3,'L':-4,'N':0,'Q':-2,'P':-2,'S':0,'R':-2,'T':-2,'W':-2,'V':-3,'Y':-3},'F':{'G':-3,'F':6,'I':0,'H':-1,'K':-3,'M':0,'L':0,'N':-3,'Q':-3,'P':-4,'S':-2,'R':-3,'T':-2,'W':1,'V':-1,'Y':3},'I':{'I':4,'K':-3,'M':1,'L':2,'N':-3,'Q':-3,'P':-3,'S':-2,'R':-3,'T':-1,'W':-3,'V':3,'Y':-1},'H':{'I':-3,'H':8,'K':-1,'M':-2,'L':-3,'N':1,'Q':0,'P':-2,'S':-1,'R':0,'T':-2,'W':-2,'V':-3,'Y':2},'K':{'K':5,'M':-1,'L':-2,'N':0,'Q':1,'P':-1,'S':0,'R':2,'T':-1,'W':-3,'V':-2,'Y':-2},'M':{'M':5,'N':-2,'Q':0,'P':-2,'S':-1,'R':-1,'T':-1,'W':-1,'V':1,'Y':-1},'L':{'M':2,'L':4,'N':-3,'Q':-2,'P':-3,'S':-2,'R':-2,'T':-1,'W':-2,'V':1,'Y':-1},'N':{'N':6,'Q':0,'P':-2,'S':1,'R':0,'T':0,'W':-4,'V':-3,'Y':-2},'Q':{'Q':5,'S':0,'R':1,'T':-1,'W':-2,'V':-2,'Y':-1},'P':{'Q':-1,'P':7,'S':-1,'R':-2,'T':-1,'W':-4,'V':-2,'Y':-3},'S':{'Y':-2,'S':4,'T':1,'W':-3,'V':-2},'R':{'S':-1,'R':5,'T':-1,'W':-3,'V':-3,'Y':-2},'T':{'Y':-2,'T':5,'W':-2,'V':0},'W':{'Y':2,'W':11},'V':{'Y':-1,'W':-3,'V':4},'Y':{'Y':7}}
		
		lowest  = INFINITY
		
		for d in data.values():
			_min = min(d.values())
			
			if(_min < lowest):
				lowest = _min
		
		for char1 in alphabet:
			self.scoring_matrix[char1] = {}
			for char2 in alphabet:
				if(char1 == "-" or char2 == "-"):
					self.scoring_matrix[char1][char2] = gap_penalty
				else:
					try:												# Save the full matrix instead of the lower triangle.
						score = data[min(char1,char2)][max(char1,char2)]
						
						# Normalization:
						score = (1.0*score/-lowest)+1
						
						self.scoring_matrix[char1][char2] = score
						self.scoring_matrix[char2][char1] = score
					except:												# (rare) Amino acid which are in alphabet but not in BLOSUM62:
						if(not self.scoring_matrix.has_key(char1)):
							self.scoring_matrix[char2] = {}
						
						if(not self.scoring_matrix.has_key(char2)):
							self.scoring_matrix[char2] = {}
						
						self.scoring_matrix[char1][char2] = mismatch_penalty
						self.scoring_matrix[char2][char1] = mismatch_penalty
		
		self.gap_penalty = gap_penalty
		
		self.disabled = False
	
	def find_score(self,amino_acid1,amino_acid2):
		"""Finds the substitution score for the substitution of two amino acids.
		
		---
		@param amino_acid1: (char) An amino acid
		@param amino_acid2: (char) Another amino acid, where one of the amino acids is substituted by the other
		
		@return: The score in the substition matrix for replacing amino_acid1 for amino_acid2
		@rtype: float or integer
		"""
		
		if(self.disabled):
			raise RuntimeError("The ScoringMatrix is not properly set, define the scores first?")
		else:
			try:
				return self.scoring_matrix[amino_acid1][amino_acid2]
			except:
				raise RuntimeError("Unknown amino-acid: "+amino_acid1+" => "+amino_acid2)
	
	def find_gap_penalty(self):
		"""Finds the gap penalty of the substition matrix. This value can also be obtained by self.find_score("-","-")
		
		---
		@return: The numerical value of a gap score. This value is typically negative since it functions as a penalty.
		@rtype: float
		"""
		
		if(self.disabled):
			raise RuntimeError("The ScoringMatrix is not properly set, define the scores first?")
		else:
			return self.gap_penalty
	
	def find_mismatch_penalty(self):
		"""	Finds the mismatch penalty of the substition matrix, if there is one single.
		
		---
		@return: The numerical value of a mismatch score. This value is typically negative since it functions as a penalty.
		@rtype: float
		"""
		
		if(self.disabled):
			raise RuntimeError("The ScoringMatrix is not properly set, define the scores first?")
		else:
			return self.mismatch_penalty
	
	def find_worst_score_for_amino_acid(self,amino_acid):
		"""Finds the 'worst' possible score for an amino acid where in this case worst means 'the highest possible score'; this is the worst case scenario since it has the lowest amount of dissimilarity.
		
		---
		@param amino_acid: The highest possible alignment score for the amino acid that has been replaced by mutations.
		"""
		
		choices = self.scoring_matrix[amino_acid].values()
		choices.append(self.gap_penalty)
		
		return max(choices)
	
	def find_all_detoriating_amino_acids_for_amino_acid(self,amino_acid,threshold = 0.5,alphabet = "ARNDCEQGHILKMFPSTWYV"):
		"""Finds all amino acids that - if they substitute from a specific amino acid - score higher than a certain treshold. This can be used to select substitutions that are likely to optimize.
		
		---
		@param amino_acid: The amino acid that the constraint is based on.
		@param threshold: The minimum value.
		@param alphabet: The amino acids that should be scanned for.
		
		@return: All amino acids that have a score higher or equal to the threshold, if they are substituted by 'amino_acid'.
		"""
		
		amino_acids = []
		
		for amino_acid2 in self.scoring_matrix[amino_acid]:
			if((amino_acid2 in alphabet) and (self.scoring_matrix[amino_acid][amino_acid2] >= threshold)):
				amino_acids.append(amino_acid2)
		
		return amino_acids



class DistanceAlign:
	"""Align sequences for the purpose of finding the most distant sequence, applying length-specific normalization.
	
	---
	@cvar TRAVEL_VERTICAL:    enumerate value (1)
	@cvar TRAVEL_HORIZONTAL:  enumerate value (2)
	@cvar TRAVEL_DIAGONAL:    enumerate value (3)
	@cvar TRAVEL_DIAGONAL_MM: enumerate value (4)
	"""
	
	# Constants (traceback signals):
	TRAVEL_VERTICAL    = 1
	TRAVEL_HORIZONTAL  = 2
	TRAVEL_DIAGONAL    = 3
	TRAVEL_DIAGONAL_MM = 4
	
	def __init__(self,substitution_matrix):
		"""Initiates the distance_align function by storing the substition_matrix.
		
		---
		@param substitution_matrix: A substitution matrix of type 'SubstitutionMarix' required for scoring the alignment.
		"""
		
		# Variable:
		self.substitution_matrix = substitution_matrix
	
	def align(self,sequence1,sequence2):
		"""Aligns sequences sequence1 with sequences2.
		
		@param sequence1: The target sequence
		@param sequence2: The biological sequence
		
		---
		@return: Returns the (one of the) optimal alignment(s) and it's properties:
			- item 0 (int):    The alignment score
			- item 1 (string): The alignment for sequence1
			- item 2 (string): The alignment in symbols where " " represets a gap, ":" a mismatch and "|" a match.
			- item 3 (string): The alignment for sequence1
			- item 4 (int):    Number of possible alignments that are optimal.
			- item 5 (list):   List with all positions in sequence1 that have an exact match with sequence2 in the alignment.
		@rtype: list
		"""
		
		# Variables:
		self.sequence1 = sequence1
		self.sequence2 = sequence2
		
		self.m = len(sequence1)
		self.n = len(sequence2)
		self.length_offset = max(self.m,self.n) - min(self.m,self.n)
		
		# Initialize matrices
		self.init_scoring_matrix()
		self.init_possibility_matrix()
		self.init_route_matrix()
		
		self.sequence1_matching_positions = []
		
		# Run alignment function
		self.align_filling_phase()
		return self.align_traceback_phase()
	
	def init_scoring_matrix(self):
		"""Initiates a matrix of (size m+1 by n+1), by setting every each element to 0, except for elements [i,0] and [0,j] which are set to i * gap_penalty."""
		
		self.scoring_matrix = [[]]
		
		for n in range(self.n+1):
			self.scoring_matrix[0].append(0)
		
		for m in range(1,self.m+1):
			self.scoring_matrix.append([0])
			for n in range(1,self.n+1):
				self.scoring_matrix[m].append(0)
	
	def init_possibility_matrix(self):
		"""Initiates a matrix of (size m+1 by n+1), by setting every each element to 1, except for elements [i,0] and [0,j] which are set to 0. The matrix stores the number of possible ways that lead to the optimal solution."""
		
		self.possibility_matrix = [[]]
		for n in range(self.n+1):
			self.possibility_matrix[0].append(1)
		
		for m in range(1,self.m+1):
			self.possibility_matrix.append([1])
			for n in range(1,self.n+1):
				self.possibility_matrix[m].append(0)
	
	def init_route_matrix(self):
		"""Initializes the 'route matrix; of (size m+1 by n+1), by setting every each element to -1. The matrix stores the routes for traceback."""
		
		self.route_matrix = []
		
		for m in range(0,self.m+1):
			self.route_matrix.append([])
			for n in range(0,self.n+1):
				self.route_matrix[m].append(-1)
	
	def align_filling_phase(self):
		"""Fills the matrices for the alignment in a dynamic programming manner."""
		
		m_max = self.m + 1
		n_max = self.n + 1
		
		# Fill the scoring, route- and possiblity matrix:
		for i in range(1,m_max):
			for j in range(1,n_max):
				# Find the highest score for substitution of amino acid i for j:
				add = self.substitution_matrix.find_score(self.sequence1[i-1],self.sequence2[j-1])
				match  = self.scoring_matrix[i-1][j-1] + add
				
				# Solves open ends / free end / free gap ends / gap free
				if(i == self.m):
					insert = self.scoring_matrix[i  ][j-1] + 0
				else:
					insert = self.scoring_matrix[i  ][j-1] + self.substitution_matrix.find_score("-",self.sequence2[j-1])
				
				if(j == self.n):
					delete = self.scoring_matrix[i-1][j  ] + 0
				else:
					delete = self.scoring_matrix[i-1][j  ] + self.substitution_matrix.find_score(self.sequence1[i-1],"-")
				
				self.scoring_matrix[i][j] = max(match,delete,insert)
				
				# Store traceback data:
				current_possibilities = []
				if(insert == self.scoring_matrix[i][j]):
					current_possibilities.append(self.possibility_matrix[i  ][j-1])
					self.route_matrix[i][j] = self.TRAVEL_VERTICAL
				if(delete == self.scoring_matrix[i][j]):
					current_possibilities.append(self.possibility_matrix[i-1][j  ])
					self.route_matrix[i][j] = self.TRAVEL_HORIZONTAL
				if(match  == self.scoring_matrix[i][j]):
					current_possibilities.append(self.possibility_matrix[i-1][j-1])
					if(add < 1):
						self.route_matrix[i][j] = self.TRAVEL_DIAGONAL_MM
					else:
						self.route_matrix[i][j] = self.TRAVEL_DIAGONAL
				
				# Find the number of optimal solutions for substition of amino acid i for j:
				for number_of_possiblities in current_possibilities:
					self.possibility_matrix[i][j] += number_of_possiblities
	
	def align_traceback_phase(self):
		"""Runs the traceback phase of the dynaminc programing algorithm; finds one of the optimal alignment of sequence1 and sequence2.
		
		---
		@return: Returns the (one of the) optimal alignment(s) and it's properties:
			- item 0 (int):    The alignment score
			- item 1 (string): The alignment for sequence1
			- item 2 (string): The alignment in symbols where " " represets a gap, ":" a mismatch and "|" a match.
			- item 3 (string): The alignment for sequence1
			- item 4 (int):    Number of possible alignments that are optimal.
			- item 5 (list):   List with all positions in sequence1 that have an exact match with sequence2 in the alignment.
		@rtype: list
		"""
		
		# Define variables:
		align1 = ""
		symbol = ""
		align2 = ""
		
		i = self.m
		j = self.n
		
		# Trace back (reverse):
		while i > 0 and j > 0:
			if(self.route_matrix[i][j] == self.TRAVEL_DIAGONAL):
				i -= 1
				j -= 1
				
				align1 += self.sequence1[i]
				symbol += "|"
				align2 += self.sequence2[j]
				
				self.sequence1_matching_positions.append(i)
			
			if(self.route_matrix[i][j] == self.TRAVEL_DIAGONAL_MM):
				i -= 1
				j -= 1
				
				align1 += self.sequence1[i]
				symbol += ":"
				align2 += self.sequence2[j]
			
			elif(self.route_matrix[i][j] == self.TRAVEL_HORIZONTAL):
				i -= 1
				
				align1 += self.sequence1[i]
				symbol += " "
				align2 += "-"
			
			elif(self.route_matrix[i][j] == self.TRAVEL_VERTICAL):
				j -= 1
				
				align1 += "-"
				symbol += " "
				align2 += self.sequence2[j]
		
		# Finish gap ends:
		while i > 0:
			i -= 1
			
			align1 += self.sequence1[i]
			symbol += " "
			align2 += "-"
		
		while j > 0:
			j -= 1
			
			align1 += "-"
			symbol += " "
			align2 += self.sequence2[j]
		
		return [self.scoring_matrix[self.m][self.n],align1[::-1],symbol[::-1],align2[::-1],self.possibility_matrix[self.m][self.n],self.sequence1_matching_positions]



class BiologicalSequence:
	"""Stores and takes care of a biological sequence and it's alignments."""
	
	def __init__(self,sequence):
		"""Initiates a BiologicalSequence by resetting all variables.
		
		---
		@param sequence: The biological sequence.
		"""
		
		self.sequence = sequence
		self.alignment = False
		
		self.reset_missed_permutations()
	
	def has_alignment(self):
		"""Checks if there has been an alignment done with this sequence and the target sequence.
		
		---
		@return: Whether there are alignment results with the target sequence stored in this object.
		@rtype: boolean
		"""
		
		if(self.alignment):
			return True
		else:
			return False
	
	def get_sequence(self):
		"""Returns the biological sequence.
		
		---
		@return: The biological sequence.
		@rtype: string
		"""
		
		return self.sequence
	
	def get_missed_mutations(self):
		"""Gives the number of missed mutations of the target sequence since the last alignment of this biological sequence with the target sequence. In the case of BLOSUM scoring not the actual number of mutations is used but the worst-case scoring scenario."""
		
		return self.missed_mutations
	
	def increase_missed_mutations(self,increase=1):
		"""Increases the number of missed mutations if this biological sequence was not in the SCOPE. Sometimes the optimization algorithms decides to apply more than one mutation in one iteration. For this reason the amount of increment can be chosen.
		
		---
		@param increase: The number of mutations to increase. Can be a float for BLOSUM62 scores.
		"""
		
		self.missed_mutations += increase
	
	def reset_missed_permutations(self):
		"""Resets the number of missed mutions. This should typically be done when the sequence was in the SCOPE and the sequence optimized the fitness function."""
		
		self.missed_mutations = 0
	
	def get_score(self):
		"""Gives the score of the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[0]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[0]")
	
	def get_alignment1(self):
		"""Gives the 'alignment1' (which is the representation of the target sequence) of the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[1]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[1]")
	
	def get_symbol(self):
		"""Gives the 'symbol' (which is the representation of similar and dissimilar amino acids) of the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[2]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[2]")
	
	def get_alignment2(self):
		"""Gives the 'alignment2' (which is the representation of the biological sequence) of the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[3]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[3]")
	
	def get_number_of_optimal_solutions(self):
		"""Gives the number of optimal solutions (log-scale) of the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[4]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[4]")
	
	def get_matching_positions_of_target_sequence(self):
		"""Gives a list with all positions in sequence1 that have an exact match with sequence2 in the last alignment of this biological sequence with the target sequence."""
		
		if(self.has_alignment()):
			return self.alignment[5]
		else:
			raise RuntimeError("Call to undefined property: BiologicSequence.alignment[5]")
	
	def set_alignment(self,alignment_results):
		"""Saves the results from an alignment with the target sequence with this biological sequence.
		
		---
		@param alignment_results: Results from the DistanceAlign::align() function.
		"""
		
		self.alignment = alignment_results



class MostDistantSequenceGenerator:
	"""Tries to find a trajectorty of sequences that have a decreased overall similarity with all biological sequences. The dissimilarity is measured by minimizing the DistanceAlign function using a specified SubstitutionMatrix. During runtime the results are directly saved to a file."""
	
	def __init__(self,job_name,start_sequence,library_sequences,substitution_matrix):
		"""Initializes the distant_sequence_generator:
			- Set all variables
			- Load biological and query sequences from file
		
		---
		@param substitution_matrix: The substitution matrix that used by the alignment function to score sequence similarity.
		"""
		
		# Define general variables:
		self.job_name = job_name
		self.filename = prefix+"most_distant_sequence__"+job_name+".txt"
		
		self.substitution_matrix = substitution_matrix
		
		self.target_sequence = start_sequence
		self.library_sequences = []
		
		self.roottime = time.time()
		
		# Define scoring variables:
		self._scoring__aligner = DistanceAlign(self.substitution_matrix)
		
		self._scoring__highest_score = False							# The highest score
		self._scoring__maximum_possibilities = False					# The highest number of possible solutions for that alignment
		
		self._scoring__higest_scoring_sequences = False					# List of highly scoring sequences. >> Length of this one makes "self._scoring__highest_score__number_of_sequences" useless
		
		# Define iteration variables:
		self._iteration__iteration = 1
		self._iteration__first_iteration = True							# Can this one be skipped?
		self._iteration__updated_previous_time = True
		self._iteration__iterions_away_from_last_succes = 0
		self._iteration__mutations_previous_time = 0
		self._iteration__mutations_previous_time__max_score = 0
		
		# Initiator functions:
		self.set_library_sequences(library_sequences)
		self.read_previous_results()
	
	def read_previous_results(self):
		"""Read the previous results, if any exist.
		
		---
		@return: If there was previous information found in "self.filename"
		@rtype: boolean
		"""
		
		try:
			found_previous_results = False
			fh = open(self.filename,"a+")
			
			if(fh):
				lines = fh.read().strip().split("\n")
				params = lines[-1].split("\t")
				
				if(len(params) > 2 and params[2] != "" and params[2].lower() != "sequence"):
					self.target_sequence = params[2]
					self._scoring__highest_score = float(params[1])
					self._scoring__maximum_possibilities = float(params[3])
					found_previous_results = True
			
			if(not found_previous_results and len(params) == 1):
				fh.write("time\tscore\tsequence\tpossiblities\n")
				found_previous_results = False
			
			fh.close()
			return found_previous_results
		
		except:
			raise RuntimeError("Couldn't access file: "+self.filename)
	
	def set_library_sequences(self,library_sequences):
		"""Store all the biological (library) sequences in the class in as an "BiologicalSequence" object.
		
		---
		@param library_sequences: The library of biological sequences.
		@type list of strings
		"""
		
		del(self.library_sequences)
		self.library_sequences = []
		
		for sequence in library_sequences:
			seq = BiologicalSequence(sequence)
			
			self.library_sequences.append(seq)
	
	def mutate_sequence(self,alphabet = "ARNDCEQGHILKMFPSTWYV"):
		"""Mutates the target sequence to a temporary new target sequence. Substition is done with the chars in 'alphabet' with a uniform probability. A selection is applied to reduce the number of non-fit mutations.
		
		---
		@param alphabet: The alphabet of which the mutations can be chosen.
		"""
		
		tmp_sequence = self.target_sequence
		
		# Define the positions in the target sequence that have similarity with biological sequences
		items = []
		for key in self._scoring__higest_scoring_sequences:
			for item in self.library_sequences[key].get_matching_positions_of_target_sequence():
				items.append(item)
		
		if(len(items) == 0):
			return False												# No mutations that improve the score are possible - reached the most distant sequence
		
		# Check how many unique positions this are and pick the square root of this number
		possible_mutations = len(sorted(set(items)))
		max_allowed_mutations = int(round(math.sqrt(possible_mutations)))
		
		# Pick randomly a number between 1 and the minumum of (square root of possiblities & arbitrary number=5) - (The reason of 5 is that no too large jumps are made: results in asmoother trajectory)
		self._iteration__mutations_previous_time = random.randint(1,min(max_allowed_mutations,5))
		self._iteration__mutations_previous_time__max_score = 0
		
		print " * Mutating "+str(self._iteration__mutations_previous_time)+" times."
		
		# Calculate which substitions should be done:
		i = 0
		while i < self._iteration__mutations_previous_time:
			position = random.sample(items,1)[0]						# Shuffle to preserve probability
			char1 = self.target_sequence[position]
			items = [x for x in items if x != position]					# Remove duplicates
			
			current_alphabet = alphabet.replace(char1,"")
			
			for key in self._scoring__higest_scoring_sequences:
				char2 = self.mutate_sequence__find_smart_position(key,position)
				
				if(char2):
					# Seleection criterion for amino acid substition
					for char in self.substitution_matrix.find_all_detoriating_amino_acids_for_amino_acid(char2,self.substitution_matrix.find_score(char1,char2)):
						current_alphabet = current_alphabet.replace(char,"")
			
			if(current_alphabet != ""):
				# Perform the actual mutation:
				perm = random.randint(0,len(current_alphabet)-1)
				char2 = current_alphabet[perm]
				
				tmp_sequence = tmp_sequence[:position] + char2 + tmp_sequence[position+1:]
				
				self._iteration__mutations_previous_time__max_score += self.substitution_matrix.find_worst_score_for_amino_acid(char2)
				i += 1
		
		self.tmp_sequence = tmp_sequence
		return True
	
	def mutate_sequence__find_smart_position(self,sequence_id,position):
		"""Finds the char in the biological sequence that corresponds in the to the target_sequence at position 'position', if any exists.
		
		---
		@param sequence_id: The id of the biological sequence in the library
		@param position: The position in the target_sequence (range: 0,...,target_sequence_length-1)
		
		@return: The char in the biological sequence that corresponds in the alignment
		"""
		
		alignment1 = self.library_sequences[sequence_id].get_alignment1()
		alignment2 = self.library_sequences[sequence_id].get_alignment2()
		symbol = self.library_sequences[sequence_id].get_symbol()
		
		if(position < 0):
			return False
		elif(position >= len(alignment1.replace("-",""))):
			return False
		else:
			seq1 = alignment1.replace("-","")
			
			diff = 0
			i = 0
			
			while(i <= position ):
				char = alignment1[i+diff]
				if(char == "-"):
					diff += 1
				else:
					i += 1
			
			pos = position + diff
			
			if(alignment1[pos] == "-"):
				return False
			
			if(pos < len(alignment2)):
				char1 = alignment2[pos]
				char2 = symbol[pos]
				if(char == "-" or char2 == ":"):
					return False
				else:
					return char
			else:
				return False
	
	def find_scope_and_calculate_fitness(self,tmp_sequence):
		"""Finds the scope (sequences of which the distance must be calculated) and calculates if the novel sequence indeed improves the fitness function.
		
		---
		@param tmp_sequence: The modified version of the target sequence which should be testet for score improvement / optimization.
		
		@return: If this iteration was successfull (if fitness was optimized).
		"""
		
		scores = {}
		gain = 0
		
		# Find scope and update
		scope_sequences = {}
		
		for key in range(len(self.library_sequences)):					# Define the SCOPE here:
			
			# The number of possible missed scores since the sequences previous alignment: previous iterations + current iteration
			missed_permutations = (self.library_sequences[key].get_missed_mutations()+self._iteration__mutations_previous_time__max_score)
			
			# Define the SCOPE criterion - optimzation stuff
			if(self.library_sequences[key].has_alignment()):
				cur_diff = (self._scoring__highest_score - self.library_sequences[key].get_score())
			
			if(not self.library_sequences[key].has_alignment() or (cur_diff <= missed_permutations)):# If is in SCOPE; either by distance or by first iteration
				scope_sequences[key] = {}
				alignment = self._scoring__aligner.align(self.tmp_sequence,self.library_sequences[key].get_sequence())
				
				print "# "+str(key+1)+" / "+str(len(self.library_sequences))+" = "+str(alignment[0])+" (steps ago since last alignment: "+str(missed_permutations)+")"
				
				# Find difference between previous score, if there is a previous alignment
				if(self.library_sequences[key].has_alignment()):
					gain += self.library_sequences[key].get_score() - alignment[0]
				
				# Remember stats - updating should be done when we're sure this indeed is an improvent of the target-sequence
				scope_sequences[key] = alignment
				scope_sequences[key][4] = math.log(scope_sequences[key][4])
				
				# Add it to the processed list
				if(not scores.has_key(alignment[0])):
					scores[alignment[0]] = []
				
				# Save the sequences ordered by score
				scores[alignment[0]].append(key)
		
		# Highest score:
		scores_sorted = sorted(scores.keys())[::-1]
		highest_score = scores_sorted[0]
		highest_score__number_of_sequences = len(scores[highest_score])
		
		# Highest possibilities:
		current_number_of_solutions = -INFINITY
		for sequence_id in scores[highest_score]:
			if(scope_sequences[sequence_id][4] > current_number_of_solutions):
				current_number_of_solutions = scope_sequences[sequence_id][4]
		
		current_number_of_solutions = scope_sequences[scores[highest_score][0]][4]
		
		# Define selection criterion(s):
		criterion0 = (type(self._scoring__highest_score) == type(False)) or (type(self._scoring__higest_scoring_sequences) == type(False)) or (type(self._scoring__maximum_possibilities) == type(False))
		criterion1 = (highest_score < self._scoring__highest_score)
		criterion2 = (highest_score == self._scoring__highest_score and (type(self._scoring__higest_scoring_sequences) != type(False)) and (
						gain > 0 or
						current_number_of_solutions < self._scoring__maximum_possibilities or
						highest_score__number_of_sequences < len(self._scoring__higest_scoring_sequences)))
		
		if criterion0 or criterion1 or criterion2:
			if(criterion0):
				print " + Sequence changed because no previous sequence was reported"
			elif(criterion1):
				print " + Sequence changed because score changed to: "+str(highest_score)
			elif(criterion2):
				if(gain > 0):
					print " + Sequence changed because a gain of: "+str(gain)
				elif(current_number_of_solutions < self._scoring__maximum_possibilities):
					print " + Sequence changed because possiblities changed from "+str(self._scoring__maximum_possibilities)+" to "+str(current_number_of_solutions)
				elif(highest_score__number_of_sequences < len(self._scoring__higest_scoring_sequences)):
					print " + Sequence changed because the number of highly scoring sequences dropped from "+str(len(self._scoring__higest_scoring_sequences))+" to "+str(highest_score__number_of_sequences)
			
			self.set_sequence(self.tmp_sequence,highest_score,scores[highest_score],current_number_of_solutions)
			
			# If object is in the scope, reset all settings, otherwise increase the number of missed mutations
			for key in range(len(self.library_sequences)):
				if(key in scope_sequences.keys()):
					self.library_sequences[key].set_alignment(scope_sequences[key])
					self.library_sequences[key].reset_missed_permutations()
				else:
					self.library_sequences[key].increase_missed_mutations(self._iteration__mutations_previous_time__max_score)
			
			return True
		else:
			print " - Score was "+str(highest_score)+" instead of < "+str(self._scoring__highest_score)
			
			# If something went wrong with the scope: reinit.
			if(self._iteration__iterions_away_from_last_succes >= 25):
				self._scoring__highest_score = False
				print "Reset..!"
			
			return False
	
	def run(self):
		"""Starts the loop and continuous as long as the time allows."""
		
		go = True
		self.tmp_sequence = self.target_sequence
		
		while(go and (abs(time.time()-self.roottime) <= 60*60*24*6)):# If it still works, or time exceeds 4 days?
			go = self.iteration()
	
	def iteration(self):
		"""Runs one iteration of the optimization algorithm; one step to valide the most_distant_sequence."""
		
		print "["+str(self._iteration__iteration)+"] # Iteration"
		print "TARGET SEQUENCE: "+self.tmp_sequence
		
		# 1) Do the comparison
		self._iteration__updated_previous_time = self.find_scope_and_calculate_fitness(self.tmp_sequence)
		
		# 2) Permute for next iteration
		permuted = self.mutate_sequence()
		
		# Complete the iteration
		self._iteration__iteration += 1
		
		return permuted
	
	def set_sequence(self,sequence,highest_score,higest_scoring_sequences,current_number_of_solutions):
		"""Saves the novel target sequence if the fitness indeed has increased.
		
		---
		@param sequence: Novel target sequence
		@param highest_score: Score of the highest sequence
		@param higest_scoring_sequences: All the highest scoring sequences
		@param current_number_of_solutions: Number of possible alignments
		"""
		
		# Store information:
		self.target_sequence = sequence
		self._scoring__highest_score = highest_score
		self._scoring__maximum_possibilities = current_number_of_solutions
		self._scoring__higest_scoring_sequences = higest_scoring_sequences
		
		# Save results:
		fh = open(self.filename,"a")
		fh.write(str(datetime.datetime.now()).split(".")[0]+"\t"+str(self._scoring__highest_score)+"\t"+str(self.target_sequence)+"\t"+str(self._scoring__maximum_possibilities)+"\n")
		fh.close()



class IO:
	"""Takes care of the in- and output of the program
	
	Program arguments:
		- What is your sequence libfile?
			* "library.txt"
		- What is your desired start-sequence?
			* Poly-W
			* First sequence of library
			* Arbitrary
		- What scorig method do you want to use?
			* Equal scoring for every AA substition
			* Proportional to frequences
		"""
	def __init__(self):
		"""Initiate all variables for in- and output."""
		
		self.library = []
		self.library_file = prefix+"library.txt"
		self.start_sequence = 1
		self.scoring_method = 1
		self.job_name = False
		
		self.parse_arguments()
		self.read_sequence_library()
		
		# Validate if the arguments were passed through correctly:
		if(not self.job_name):
			self.usage()
		
		# Create substitution matrix accoriding to the parameters:
		substitution_matrix = SubstitutionMatrix()
		if(self.scoring_method == 1):
			substitution_matrix.set_equal_scoring_probabilities()
		else:
			substitution_matrix.set_blosum62()
		
		# Draw initization settings:
		print "# Start sequence: "+str(self.start_sequence)
		print "# Scoring method: "+str(self.scoring_method)
		print "# Library file:   "+str(self.library_file)
		print
		
		# Run the optimization algorithm with the correct parameters:
		most_distant_sequence_generator = MostDistantSequenceGenerator(self.job_name,self.start_sequence,self.library,substitution_matrix)
		most_distant_sequence_generator.run()
	
	def parse_arguments(self):
		"""Parses the arguments. Details about the arguments are given in IO::usage()"""
		
		l = len(sys.argv[1:])
		i = 1
		
		while(i <= l):
			arg = sys.argv[i]
			if(arg == "-h" or arg == "--help"):							# -h --help
				self.usage()
			elif(arg == "-v" or arg == "--version"):					# -v --version
				self.version()
			elif(arg == "-i"):											# -i			Initiator sequence
				i += 1
				if(i <= l):
					cmd = sys.argv[i]
					if(cmd == "1"):
						self.start_sequence = 1
					elif(cmd == "2"):
						self.start_sequence = ""
						for j in range(300):
							self.start_sequence += "W"
					else:
						self.start_sequence = cmd
				else:
					self.usage()
			elif(arg == "-s"):											# -s			Scoring function
				i += 1
				if(i <= l):
					cmd = sys.argv[i]
					if(cmd == "1"):
						self.scoring_method = 1
					elif(cmd == "2"):
						self.scoring_method = 2
					else:
						usage()
					
				else:
					self.usage()
			elif(arg == "-f" or arg == "-l"):							# -f -l			File / Library
				i += 1
				if(i <= l):
					cmd = sys.argv[i]
					self.library_file = cmd
				else:
					self.usage()
			else:														#				Jobname
				if(i == 1):
					self.job_name = sys.argv[1]
				else:
					usage()
			
			i += 1
	
	def read_sequence_library(self):
		"""Reads all sequences from the file "self.library". The file is of the following format:
		"SEQUENCE1\\nSEQUENCE2\\nSEQUENCE3"
		"""
		
		try:
			fh = open(self.library_file,"r")
			self.library = fh.read().strip().split("\n")
			fh.close()
			
			if(self.start_sequence == 1):
				self.start_sequence = self.library[0]
		
		except:
			raise RuntimeError("No access to file: "+self.library_file+" or file is invalid")
	
	def usage(self):
		"""Prints the usage."""
		
		print "Usage: "+sys.argv[0]+" job_name [arguments]"
		print
		print 
		print "The following parameters can be used:"
		print "  -f LIBRARY_FILE   File containing biological sequences; newline delimited"
		print "                    'library.txt' (Default)"
		print "  -i 1/2/SEQ        Initiator sequence:"
		print "                    1   First sequence in the library (Default)"
		print "                    2   Poly W (len=300)"
		print "                    SEQ A custom sequence"
		print "  -s 1/2            Scoring function of amino acids:"
		print "                    1   Equal / uniform (Default)"
		print "                    2   BLOSUM62"
		print "  -v                Shows the version and license"
		print
		print
		print "Please send your bug-reports to:"
		print "   <https://code.google.com/p/yh-most-distant-sequence/issues/list>"
		sys.exit()
	
	def version(self):
		"""Prints the scipts' version and additional LICENSE information."""
		
		print "[Version]"
		print "  yh-most-distant-sequence: "+__version__
		print
		print "[License] :: "+__license__
		print "  Copyright (c) 2013 ing. Youri Hoogstrate"
		print
		print "  Permission is hereby granted, free of charge, to any person"
		print "  obtaining a copy of this software and associated documentation"
		print "  files (the \"Software\"), to deal in the Software without"
		print "  restriction, including without limitation the rights to use,"
		print "  copy, modify, merge, publish, distribute, sublicense, and/or sell"
		print "  copies of the Software, and to permit persons to whom the"
		print "  Software is furnished to do so, subject to the following"
		print "  conditions:"
		print
		print "  The above copyright notice and this permission notice shall be"
		print "  included in all copies or substantial portions of the Software."
		print
		print "  THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND,"
		print "  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES"
		print "  OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND"
		print "  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT"
		print "  HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,"
		print "  WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING"
		print "  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR"
		print "  OTHER DEALINGS IN THE SOFTWARE."
		print
		print "[Bugs]"
		print "  Please send  your bug-reports, feedback and questions to:"
		print "    <https://code.google.com/p/yh-most-distant-sequence/issues/list>"
		
		sys.exit()

prefix = '/home/yhoogstrate/most_distant_sequence/'
io = IO()

print "Done"
